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1 INTRODUCTION 



Recently, ID' Alessio etail d2005h ; ICalvet et al.1 J2005T) identified 
three young stars (CoKu Tau/4, GM Aur and DM Tau) where the 
circumstellar disc appears to have an inner hole. These observa- 
tions were performed with the Spitzer Space Telescope. CoKuTau/4 
shows no evidence of continued accretion, while GM Aur and DM 
Tau are still accreting. GM Aur even shows some evidence for 
an optically thin disc in the innermost portions of the hole. These 
systems are often referred to as 'transition discs .' It has been sug- 
1 gested that these holes wer e cleared by a planet JOuillen et alj2004l ; 
IVarniere et al]|2004l2006h . 

In this paper we shall discuss the 3d structure of gaps cre- 
ated in discs by embedded planets, and the possible observational 
consequences. This is important, since Spitzer cannot image these 
systems directly; we are obliged to analyse spectral energy distribu- 
tions (SEDs). A key problem in interpreting Spitzer measurements 
is our lack of knowledge about the dust and gas dynamics close to 
a planet-induced disc gap. In particular, the disc edges appear to 
have h/r ~ 0.1, which is too large for simple (two dimensional) 
thin disc models to be wholly ap propriate. 

Despite these limitations, iRice et al. I d2006h have already 
demonstrated (using 2d calculations) the interesting possibil- 
ity of a 'dust filter' acting close to the edge. This could be 
a way of permitting GM Aur and DM Tau to continue to ac- 
crete (cf iLubow and D' Angeloll2006l. below), while still retaining 
their disc holes. Several authors ([ Paardekooper and Mellem ai2004l ; 
Madd ison et alj|2007l ; iFouchet et alj|2007l) have noted that a gap 
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ABSTRACT 

Giant planets embedded in circumstellar discs are expected to open gaps in these discs. We 
examine the vertical structure of the gap edges. We find that the planet excites spiral arms with 
significant (Mach number of a half) vertical motion of the gas, and discuss the implications 
of these motions. In particular, the spiral arms will induce strong vertical stirring of the dust, 
making the edge appeared 'puffed up' relative to the bulk of the disc. Infra-red observations 
(sensitive to dust) would be dominated by the light from the thick inner edge of the disc. 
Sub-millimetre observations (sensitive to gas velocities) would appear to be hot in 'turbulent' 
motions (actually the ordered motion caused by the passage of the spiral arms), but cold in 
chemistry. Resolved sub-millimetre maps of circumstellar discs might even be able to detect 
the spiral arms directly. 

Key words: planetary systems : protoplanetary disks 



might open in the dust, before the gas disc shows similar behaviour. 
This is due to the dust dis c having zero viscosity . Vertical settling 
of the dust within the disc jD' Alessio et al.ll2006h will also be rele- 
vant. If any cross-gap flow is primarily from the upper layers of the 
disc, then it will naturally be depleted in dust. 

Most previous numerical w ork on protoplanet formin g discs 
has been performed in 2D (see de Val-Borro et al. 2006, for a 
selection of codes which have been used to study the problem), 
due to the large computational cost of 3d calculations. Most pre- 
vious work in 3d has conc entrated on migration (e.g. I Kiev et all 
l200ll : lD'Angelo et alj|2002l). a nd the flow in the circumplanetary 
region faahr and Kiev! l2006h . iBate et all d2003l) studied migra- 
tion rates of planets in 3d discs, finding tha t 3d consideration s 
slowed Type I migration rat e s (as p redicted bv lTanaka et al.l 2002). 
Paardekooper and Me llemal fc006l) studied the effect of including 
radiative transfer in their 3d code, concluding that the pressure gra- 
dients induced could reverse Type I migration. These studies have 
illustrated that 2d calculations are poor approximations to the true 
flow in the vicinity of the planet. Of course, for global calculations 
of the effect of lower mass planets on a disc, the 2d approximation 
is more than adequate (and we shall see an aspect of this in our re- 
sults below). In this paper, we shall demonstrate important effects 
in the gap edges which are only seen in 3d calculations. 

iBolev et al.l fc005l) ; lBolev and Duri sen ( 2006) studied the ver- 
tical structure of shocks in self-gravitating discs, without an em- 
bedded protoplanet They particularly examined the behaviour of 
the gas before and after the shock jump. Depending on the nature 
of the equation of state and the degree to which self gravity was 
important, they found that the disc might compress or expand ver- 
tically after passing a shock. They find that a breaking wave could 
be excited on the surface layers of the disc, which has the potential 
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to mix the disc radially. In this paper, we discuss the effect of the 
spiral arms raised by an embedded planet. 

In this paper, we discuss our numerical method in Section[2] 
We present our results in Section|3] and discuss their consequences 
in Section|4] Our conclusions are presented in Section|5] 



2 NUMERICAL SETUP 

In this section, we shall describe the code we used in our numerical 
experiments. Full 3d models of circumstellar discs are extremely 
challenging, requiring over an order of magnitude more computing 

power than 2d models of comparable resolution. 

Our code is based on the FLASH code of iFrvxell et al.l 

(2000), an adaptive mesh refinement (AMR) code based around 
a Piecewise-Parabolic Method (PPM) hydrodynamics solverQ We 
made several modifications to the code, to adapt it for these ex- 
periments. First, we adjusted the transport step to enforce con - 
servation of angular momentum about the z-axis dKlevlll998h . 
The detailed implentation of this is that used by the polar FLASH 
code in the comparison of hydrodynamics codes presented by 
Ide Val-Borro et alj {2006). We have verified that the z component 
of the angular momentum is conserved to machine precision. Our 
discs were isothermal in the azimuthal and vertical directions, and 
had a temperature gradient in the radial direction. The pressure at 
any point in the disc is calculated from this temperature profile and 
the local density, providing the input states to the Riemann solver. 
In solving the Riemann problem, the solver assumes that the gas has 
7 = 1.1. The radial temperature gradient is set to produce a flaring 
disc (h/r oc r a for some a). We added a simple NBODY solver, to 
model the presence of a planetary system. In this paper, we only 
use two bodies (a star and a planet), and they do not feel the gravi- 
tational effect of the disc. Consequently, the orbit of the planet did 
not evolve in runs presented herein. In addition, we added physical 
viscosity to the code. We implemented all components of the vis- 
cous stress tensor, but the calculations shown here only use the r-(j> 
component, since it is not clear that the true effects of MH D turbu- 
lence can be reproduced by a physical viscosity (see I Winters et al.l 
120031) . Buffer zones in the radial direction damp waves, preventing 
them reaching the inner or ou ter boundaries. The buffe r zones are 
implemented in the manner of Ide Val-Borro et alj J2006T) . 

Our computational grid ranges between r = 2 x 10 13 cm and 
2 x 10 14 cm, between z — and 3 x 10 13 cm vertically and over 
a full circle in azimuth. The buffer zones reduce the effective radial 
extent to 3.2 x 10 13 cm < r < 1.7 x 10 14 cm. We used a uni- 
form mesh, with (n r ,n^,n z ) = (128,384,64). This makes the 
grid cells almost square around the planet in the r-cj> plane, but 
with higher resolution in the z direction. The increased 2 resolu- 
tion is necessary to avoid under-resolving the vertical structure of 
the inner disc. We do not use the AMR capabilities of FLASH, since 
we are not certain that the jumps in refinement do not produce any 
artifacts. The edge of the gap will always lie roughly parallel to the 
fine-coarse interfaces, opening the possibility that small errors in 
the interpolation between levels will b uild up over time (see also 
the comments o f lCiecielag etafll 20071) . Our planetary system con- 
sisted of a planet in a fixed circular orbit around a 1 Mq star. We 
performed runs with planetary masses of 1 and 2 Mj. The orbital 
semi-major axis was 7.5 x 10 13 cm = 5 au. 

The initial density profile of the disc was constructed so that 



the surface density was constant at Eo = 10 3 g cm -2 . There is 
great uncertainty in the appropriate density profile of protoplane- 
tary discs, with estimates ranging between the constant value we 
use here, and the E oc r -3 / 2 form of the Minimum Mass Solar 
Nebula. The disc was mildly flaring, with 
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with rp = 1.5 x 10 cm. We made use of equation 4 of 
iTanaka et ai] d2002l) to adjust the initial orbital velocity of the gas 
to take account of the radial pressure gradient. 

We neglect self-gravity in our calculations. Together with our 
assumption of a vertically isothermal disc, th is prevents u s repro - 
ducing much of the beha viour observed by iBolev et all d2005l) : 
iBolev and Duriser] d2006t) . However, since our discs contain a 
planet, self-gravity is not necessary to generate spiral arms. The 
assumption of a vertical isothermali ty also supresse s the 'wave 
channelling' behaviour described bv lBate et al.1 d2002h . where spi- 
ral waves refract towards the surface of a vertically stratified disc. 
Nevertheless, these calculations provide a useful starting point for 
analysing the vertical structure of gap edges. 



3 RESULTS 

We shall present the results of our numerical experiments, grouped 
by the planet mass. For each planet mass, we performed runs with 
different values for the coefficient of viscosity, v. The first run set 
v = 0, and hence is dependent solely on the intrinsic numerical 
dissipation of the codej 2 ] Two non-zero values were tested, v = 
10 15 cm 2 s _1 and v = 10 16 cm 2 s _1 . We shall refer to these runs 
are 'zero,' 'low' and 'high' viscosity respectively. 

The high and low viscosity discs have Reynolds numbers 
at the planet's orbit of 1Z = r 2 Q,/u = 10 5 and 10 4 respec- 
tively. There are a number of different criteria for determining 
whether a planet should o pen a gap. First is the tidal criterion of 
iLin and Papaloizoul dl993l) : 



q > 3 {h/rf 



(2) 



This is obtained by comparing the size of the Hill sphere, r H = 
r(q/3) 1//3 , which is the volume over which the planet's gravity 
dominates, to the thickness of the disc. More commonly used is 
the viscous condition 

q > 4072T 1 (3) 

as discussed bv lBrvden et al. l dl999h . lCrida"et al. (2006) combined 
these two conditions into a single one: 

4 r H qTZ 

In Table Q] we show the gap-opening predictions for each criteria, 
listed by planet mass and disc viscosity. There is significant uncer- 
tainty, underlined by the fact that the defi nition of a gap is some- 
what arbitrary (see Hosseinbor et al. 20071). 

We shall study the state of the disc a fter 100 orbits. Al- 
though a stea dy state is never reached (see, e.g. lVarniere et al J2006I : 
lEdgaj|2007t) , the condition of the dis c after 100 orbits is gen - 
erally accepted as being representative dde Val-Borro et al. 2006). 



1 The source code is available at |http : //flash . uchicago ■ edu/| 



Note that there is no reason to expect the numerical dissipation to act like 
a viscous term 
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Criterion 


1 Mj 


2Mj 


Tidal 


ggg 


ggg 


Viscous 


ggn 


ggn 


Combined 


gnn 


ggn 



Table 1. Predictions of gap opening according to different criteria. The cri- 
teria are explained in the text. Each entry lists three values, for the zero, low 
and high viscosity cases respectively. A 'n' indicates that no gap should 
form, a 'g' that a gap is expected 



Each orbit requires approximately 300 h of CPU time on NCSA's 
mercury cluster. 
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3.1 Jupiter Mass Planet 

We will now discuss the results we obtained from a Jupiter mass 
(q = 10~ 3 ) planet embedded in the zero, low and high viscosity 
discs. In Figure Q] we show the surface density of the disc after 
one hundred orbits. Note that although the total surface density of 
our disc is 1000 g cm -2 , we are only computing the portion of 
the disc with z > 0. Hence the initial surface density was only 
500 g cm -2 in the computational volu me. The behaviour we se e 
is expected (see, e.g. the comparisons of Ide Val-Borro et alj[ 2006 ) . 
The zero viscosity case has lots of fine structure and a very clean 
gap. The low viscosity case has a clean gap, but most of the fine 
structure has been suppressed, while the gap in the high viscosity 
case has been suppressed (although we can still see a substantial 
depression in the surface density). 

In Figure[2]we plot the vertical Mach number, A4 Z = v z /c s 
in the disc for slices with z = r H /4 and z = r H . Unsurprisingly, 
there are strong motions close to the planet. However, we also see 
vertical motions associated with the spiral arms seen in Figure Q] 
These motions are relatively strong (up to 7V4 Z « 0.5), and reach 
well into the gap edge (although they do not propagate as far as 
the spiral arms in the density field.. The strength of these motions 
does not appear to be affected by the viscosity, with a similar struc- 
ture appearing for the zero, low and high viscosity cases. This is 
not especially surprising, since viscous effects take several orbits 
to become apparent, whereas the planet stirs each part of the disc 
every orbit. Furthermore, we only apply the r-cj> component of the 
viscous stress tensor, so the vertical motion can only be affected 
indirectly. 

Averaging the motion around an orbit, we find that the mean 
velocity is zero, indicating that the vertical motions are waves. 
Given their location, we conclude that these vertical motions are 
the 3d aspe ct of the spiral arms g enerated by the planet. Previous 
studies (see IWard and Hand 12000. and references therein) of the 
spiral arms generated by a planet have usually assumed that the 
disc is two dimensional. Here, we see that the spiral arms have a 
definite three dimensional structure. The gas contracts and expands 
vertically as it passes through the spiral arms. However, this 3d 
structure is only apparent close to be planet, as the vertical motions 
rapidly fade at larger radial distances. This can be understood in 
terms of the thin disc approximation and the planet' s Hill sphere. I f 
a disc is thin, then it may be treated as a 2d entity (Prinriej [l98lh . 
The discs used in this paper have h/r set according to EquationQ~| 
and as such are relatively thin. Accordingly, we should expect the 
waves to be 2d in nature. However, within the planet's Hill sphere 
the planet's gravity dominates. This has a strong vertical compo- 
nent, which drives the 3d structures we see in Figure|2] As these 
waves propagate beyond the planet's Hill sphere, the thin disc ap- 
proximation reasserts itself, and the vertical motions dissipate. 
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Figure 1. Surface density after 100 orbits for aq = 1Q~ li planet in the zero 
(top), low (middle) and high (bottom) viscosity cases. The scale is linear, 
and ranges between and 1000 g cm -2 . Note that since we only model 
the disc for z > 0, the unperturbed surface density is 500 g cm -2 



There are a number of possible numerical issues with our re- 
sults. Firstly, there is our choice of equation of state. We imposed 
a temperature structure on the disc, but used 7 = 1.1 in the Rie- 
mann solver. Of course, using an isothermal Riemann solver would 
not be fully self-consistent either, since there is a radial tempera- 
ture gradient. In Appendix lAl we demonstrate that reducing 7 to 
a value of 1.01 (thereby making the Riemann solver more closely 
approximate the isothermal case) does not have a significant effect 
on our results. 

The viscosity is also an issue. We chose to include only the 
r-cp component of the viscous stress tensor, which is sufficient to 
drive accretion. The true source of angular momentum transport in 
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Figure 2. Mach number of the vertical velocity on planes with z = r H /4 (left) and r H (right) for a q = 10 3 planet in the zero (top), low (middle) and high 
(bottom) viscosity cases. These are shapshots taken after 100 orbits. The scale ranges from A4 Z = — 1 (black) to 1 (white) 



astrophysical discs is unlikely to act like a simple physical viscos- 
ity, so including components not essential to accretion did not seem 
necessary. This said, it is important to verify that this does not pro- 
duce any important inconsistences. In Appendix [B] we show that 
including all components of the viscous stress tensor in our runs 
does not affect the vertical oscillations we observe. 

Finally, our cylindrical grid has relatively poor resolution 
(compared to h) in the innermost regions of the disc, so Figure [2] 
is somewhat noisy there. However, the spiral arms remain clearly 
visible, indicating that the noise is not dominating our results. Nev- 
ertheless, we performed a short resolution test, detailed in Ap- 



pendix[C] We show that doubling the resolution supresses the noise 
seen in Figure [2] with little effect on the vertical oscillations in- 
duced by the planet. 

3.2 Two Jupiter Mass Planet 

In this section, we shall discuss the results obtained when a 2 Mj 
planet orbits in the same three di scs. For this case, we fou nd that 
the buffer zone damping time (see lde Val-Borro et alj |2006) had to 
be reduced to half an orbit, to suppress wave reflection. 

We start with Figure(3] which shows the surface density in the 
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Figure 3. Surface density after 100 orbits for a q = 2 X 10 planet in 
the zero (top), low (middle) and high (bottom) viscosity cases. The scale is 
linear, and ranges between and 1000 g cm~ 2 . Note that since we only 
model the disc for z > 0, the unperturbed surface density is 500 g cm - 2 



disc after approximately 100 orbits of the planet. This can be com- 
pared to Figure Q] We see that a substantial gap has appeared in all 
cases, although there is some material visible in the gap for the high 
viscosity case. In the zero viscosity case, we can see a density en- 
hancement moving along the gap edge. When viscosity is present, 
it appears that formation of this enhancement is suppressed. As in 
the case of FigureQ] the surface density plots are much as expected. 

Figure [4] shows M z on planes with z = r H /4 and z — r H , 
and may be directly compared to Figure|2] The A4 Z arms are more 
prominent, particularly on the z = r H surface. However, the veloci- 
ties reached are not substantially greater and the radial propagation 
of the velocity field does not appear to be greater. We conclude 



that increasing the planet's mass increases the vertical domain over 
which the planet can drive vertical motions, but does not substan- 
tially increase their maximum strength or radial domain. 

The vertical domain of the vertical oscillations expands be- 
cause the driving is stronger. Doubling the mass of the planet will 
double its gravitational field, increasing the vertical distance over 
which driving is effective. However, the extent of the driving in the 
orbital plane is still limited by the size of the Hill sphere. Since 
r H oc q 1 ^ 3 , we should not expect a 2 Mj planet to excite vertical 
oscillations over a substantially larger radial domain than a 1 Mj 
planet. 

3.3 Low Mass Planet 

In addition to the two sets of runs described above, we performed 
a set of runs with a 0.1 Mj planet. As predicted, this planet did 
not open a gap, but merely caused a depression in the surface den- 
sity. Furthermore, we did not find an substantial vertical motions 
being driven, in contrast to those shown in Figures [2] and [4] This 
is consistent with the arguments about the planet's Hill sphere pre- 
sented above. The Hill sphere of a 0.1 Mj planet fits entirely within 
one vertical scale height. Since the planet's gravity only dominates 
within the Hill sphere, it is then not surprising that the vertical com- 
ponent of the planet's gravity is unable to drive strong vertical os- 
cillations. The region it affects is simply too small. 



4 DISCUSSION 

In Section[3]we presented an analysis of the vertical velocity struc- 
ture induced by a Jupiter-mass planet embedded in a circumstellar 
disc. Significant vertical motions were induced in the spiral arms 
raised by the planet, with vertical Mach numbers up to ~ 0.5. We 
shall now discuss the possible observable consequences. 

4.1 Thickness of the Gap Edge 

Vertical motions induced by the spiral arms have the potential to 
alter the appearance of the gap edge in a number of ways. Since 
the gap edge is likely to be the first feature which will be imaged 
directly, it is important to understand these. 

Firstly, if the inner disc is depleted (not in our current calcu- 
lations, but recall systems such as CoKu Tau/4, DM Tau and GM 
Aur), then stellar radiation will directly strike the gap edge, and 
heat it. This would increase the scale height, but only modestly 
since ft oc c s oc vT. The edge of the gap would also appear to 
be very bright, due to its large directly illuminated surface area. A 
more significant increase may come from the density variations in 
the spiral arms themselves. Across the spiral arm, densities vary 
by a factor of a few. However, since the density distribution in a 
vertically isothermal disc is a gaussian, the scale height only de- 
pends on the log of this. We might expect the spiral arms (if not re- 
solved) to i ncreas e the apparent scale height by w h. We note that 
I Wolf et al] 12002) made simulated observations of a similar situa- 
tion, and found that the spiral arms were difficult to see. However, 
their work appears to have been based on a 2d hydrodynamic code, 
which would suppress the strong vertical motions we found. 

Since instruments such as Spitzer observe the dust (and not 
the gas), it is important to consider the effect of the coupling be- 
tween the gas and dust. The strong vertical motions imply that gas 
streamlines move up and down about half the scale height. Dust 
coupled to gas would also move up and down by the same distance. 
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Figure 4. Mach number of the vertical velocity on planes with z = r^/A (left) and r H (right) for a q = 2 X 10 3 planet in the zero (top), low (middle) and 
high (bottom) viscosity cases. The discs are shown after 100 orbits of the planet. The scale ranges from M z = — 1 (black) to 1 (white) 



The strong vertical streaming motions imply that there are larger 
velocity gradients present than expected from a stratified or sedi- 
mented disk. We expect these streaming motions to suppress dust 
sedi mentation, as compared w i th a stan dard turbulent disc. 

IDullemond and Dominikl {2004b) developed a simple model 
for dust sedimentation, and applied it to evolving protoplanetary 
discs. In their models, dust grains larger than around 1 \im were 
sufficiently decoupled from the gas to undergo significant sedimen- 
tation, potentially affecting the expected spectral energy distribu- 
tion (SED) obtained from the discs. 

In their model. lDullemond~an d Dominik balanced the stirring 
and settling timescale of the dust. Settling was due to the vertical 
component of gravity within the disc. Balancing this against Ep- 
stein drag (which should be applicable to all dust particles within a 
protoplanetary disc), they determined a settling time of 



a E(r) 



exp 



(5) 



3 % /2lr" m fi kep \ 2h?(r 
where the grains have mass m, and cross section a. Stirring 



was provided by the viscous processes in the disc. Particles 
would diffuse upwards, under the influence of eddies in the disc, 
so the stirring timescale would occur on a timescale t stir = 
z 2 ID, where D is the e ffective diffusion coefficient. Following 
Dullemon d and Dominikl : 



Sc 



(6) 



where we have assumed that we have a standard a-disc, with 
v — ac s h, and also that D = v. Where Sc is the Schmidt number 
of the dust grains in the gas, which measures the strength of the 
coupling between the gas and the dust. In our discs, Sc only be- 
comes significantly greater than 1 around five scale heights above 
the midplane. Accordingly, we shall take Sc = 1 in the cur rent 
discussion. See section 2 of IDullemond and Dominikl d2004bl) for 
detailed exposition of dust settling. Although we used a constant 
viscosity disc, we can define an equivalent a value, finding it to be 
a w 2.9 x 10 -3 at the planet's orbital radius for the low viscosity 
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disc. The settling height is then defined by requiring t stk = £t sett . 
with £ typically one hundred. 

Our numerical experiments m ake use of a much denser gas 
disc than iDullem ond and Dominik] - see particularly their equa- 
tion 19. Their disc is depleted by a factor of roughly 30, as com- 
pared with ours. This strengthens the coupling between the dust 
and gas, i n turn increasing the time required for dust to settle to the 
midplane. iPaardekooper and Mellemal ( 120040 made use of a disc 
depleted by a similar amount, in their work showing decoupling 
between the gas and dust discs. They also followed grains of size 

I mm, which are observable to instruments such as the SMA and 
ALMA, but less easily detectable by Spitzer. We would not expect 
our disc to sediment at all. However, since we do not include self- 
gravity and impose a temperature profile, the densities in our nu- 
merical experiments can easily be rescaled. 

At later times, the disc would be depleted and we would 
expect sedimentation alo ng the lines of the theory developed by 
IDullemond and Dominikl . However, we must also consider the stir- 
ring provided by the planet. As shown in Section [3] there are 
strong vertical motions in the spiral arms within the gap edge. 
The vertical stirring caused by the spiral arms occurs on the syn- 
odic timescale, over approximately one scale height. The synodic 
timescale is given by Sl syn = Q — fi p , and varies with position in 
the disc. However, we are interested in particles close to the gap 
edge, which invariably lies close to the m — 2 Lindblad reso- 
nances. The synodic period will therefore be about twice the orbital 
period. A particle that is lifted upwared by spiral arm and then dif- 
fuses out of the arm would enter a stream line with a greater mean 
height. This could raise the effective diffusion coefficient for the 
dust grains from D w ac s h to D ss 0.5c 3 ft. Comparing to our 
equivalent a = 2.9 x 10~ 3 , we see that the stirring timescale im- 
plied by Equation [6] would fall by over two orders of magnitude. 
However, the settling time (Equation [5j is unaffected by the spiral 
arms, and so will remain the same. Because the ratio of the settling 
timescale to the stirring timescale depends on the local gas density 
(which roughly follows a gaussian distribution vertically), we ex- 
pect the depletion height to depend on the square root of the log 
of the diffusion coefficient. Since our effective diffusion coefficient 
increases by a factor of « a, we expect the depletion height to be 
raised by a few scale heights. 

In summary, three different effects lead to apparent thickening 
of edges of discs near planets 

(i) Direct stellar radiation raising the scale height and temper- 
ature in the gap edge. Size of effect: Factor of a few increases in 
temperature but scale height only depends on square root of tem- 
perature so increase in scale height is only modest 

(ii) Density variations caused by spiral arms lead to raising of 
r = 1 surface. Density contrasts in spiral arms are a factor of a few 
but variation in scale height depends on the natural log of this. Size 
of effect ~ 1 scale heights 

(iii) Increased vertical mixing of dust caused by the vertical mo- 
tions in the spiral arms. The stirring timescale is effectively in- 
creased by a~ 1 . The depletion height could vary by the square root 
of the log of this, which is a few scale heights 

Combined, these effects can make the inner edge of the gap ap- 
pear significantly thic ker than the bulk of the disc. We note that 
ID' Alessio et al.l q2005) calculated an aspect ratio of h/r w 0.1 for 
the edge of CoKuTau/4's disc, which seems quite large for a disc 

II au from a 0.5 Mq star. In addition to these effects, the vertical 
motions themselves will make some portion of the disc edge appear 
higher - perhaps by as much as half a scale height. However, since 



the vertical motions are waves, other portions of the disc edge will 
be depressed by a similar amount. Whether this extra thickening (or 
thinning) is observable will therefore depend on the viewing angle 
and orbital phase of the planet, especially in azimuthally averaged 
spectral energy distributions. 

In a fully self-consistent model, with radiative transfer and 
dust settling, the increased scale height might cause the gap to 
become shallower, as the tidal condition introduced at the start of 
Section[3]fails. It is also possible that the thick disc edge will cast a 
shadow over the outer disc, causing it to cool and b ecome geometri- 
cally thinner (see lDullemond and Dominikl 2 004a. for a discussion 
of self-shadowing in the context of Herbig Ae/Be stars). The light 
observed from such a disc would be completely dominated by the 
bright inner edge of the gap induced by the planet. 



4.2 Detection of the Velocity Field 

Almost all extra-solar planets have been discovered by the radial 
velocity (RV) method. Although its success is unquestioned, the 
RV method suffers from two strong biases. Firstly, it is most sen- 
sitive to massive planets in close orbits. Secondly, it can only be 
used on old stars (t > 1 Gyr), or surface activity will wipe out the 
small signal from a companion. We would greatly deepen our un- 
derstanding of planet formation, if we could catch planet formation 
in pro gress. New instruments, primarily ALMA, of fer this possi- 
bility. iWolf etal.1 d2002h : IWolf and D'Angeld d2005h have already 
argued that the 'accretion hotspot' of an embedded protoplanet will 
be detectable in ALMA images. We shall now discuss the prospects 
of using ALMA's high velocity resolution to detect the spiral arms 
rais ed by a protop l anet. 

iDartois et al.l d2003D examined the structure of DM Tau's outer 
(r > 50 au) accretion disc in several CO lines. They found that 
the velocity dispersion in the disc midplane was higher than that 
towards the surface. This implies that the turbulence is greater in the 
midplane, even though they found that the disc surface was warmer 
than the midplane. A natural explanation for this would be spiral 
arms stirring up the midplane. A planet is unlikely so far out in the 
disc, but gravitational instability would be quite likely. 

Even with its 0.05" resolution, the beam of ALMA will still 
be approximately 5 au across at the nearest star-forming clouds 
(located approximately 100 pc away from us). This means that any 
effects of a planet in a 5 au orbit are highly unlikely to be visible. 
Accordingly, we have rescaled our results to place the planet on 
a 10 au orbit. We computed the line of sight velocity for every 
grid cell in our computational domain, and projected this onto the 
sky. We then traced rays through the projected density structure, 
evalutating the integral 

£, os = [ pdl (7) 

J DC 

We assumed that ALMA would observe the velocity in the grid cell 
where Ei os = 1. Although this isn't necessarily the surface ALMA 
would see, it serves as a useful initial estimate. 

In Figures [5] and [6] we show two simulated position- velocity 
plots, based on a Jupiter mass planet in the low viscosity disc. The 
plots differ only in the beam track. The disc was inclined at 45° to 
the line of sight, and we average the velocities over a 5 au beam. 
The planet was on the 'far' side of the disc, in approximately the 
2 o'clock position. The beam track of Figure [5] (the 'near' side of 
the disc) does not pass over the planet, while that of Figure|6]does. 
Note that there is a clear difference in velocity structure between 
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Figure 5. Simulated ALMA observation of a Jupiter mass planet in a 10 au 
orbit. The orbital plane is inclined at 45° to the line of sight, and the planet 
is located approximately at the 2 o'clock position. The line of sight velocity 
is shown in the left panel, along with the beam track (which does not pass 
over the planet). The circle placed on beam track indicates the 5 au size 
of the ALMA beam. In the right panel, we show the position-velocity plot 
produced along the track 



putable from observations of disc chemistry. If these observations 
suggested a more reasonable temperature, then one might conclude 
that a planet was stirring up the disc, causing the mismatch between 
the 'turbulent' temperature (in reality, the spiral arm stirring) and 
the chemical temperature. Instruments such as the Submillimeter 
Array (SMA) are already making deta iled observations of circum- 
stellar discs (e.g. lOi et al ] |2004 2006), and analysis of these sys- 
tems may well find anomalies which can be explained by the pres- 
ence of a planet. 

Despite the caveats raised above, we have demonstrated the 
tantilising possibility of detecting a young planet embedded in is 
nascent disc. Such a detection would greatly extend our understand- 
ing of the planet formation process. Observational considerations 
have restricted all previous planet searches to old systems, where 
the disc has long since vanished. 




Figure 6. Simulated ALMA observation of a Jupiter mass planet in a 10 au 
orbit. The orbital plane is inclined at 45° to the line of sight, and the planet 
is located approximately at the 2 o'clock position. The line of sight velocity 
is shown in the left panel, along with the beam track (which passes over the 
planet). The circle placed on beam track indicates the 5 au size of the ALMA 
beam. In the right panel, we show the position-velocity plot calculated for 
this beam track 

the 'near' and 'far' tracks, indicating that we can easily identify 
which side of the disc we are inspecting. 

Of greater interest, is the fact that the planet is detectable 
in the right hand panel of Figure [6] as the extra contours around 
x — 10 au and y = 6 km s _1 . This is slightly outside the planet's 
actual location (after projecting onto the sky), indicating that the 
detection is of the spiral arm raised by the planet. Of course, real 
observations will be noisy, and the detection in Figure [6] is only 
in one contour. Furthermore, our approach of using the Ei os = 1 
surface to construct our velocity field is crude. Real observations 
will have to select particular molecular transitions, based on the 
temperature structure of the disc. Optically thin lines will convolve 
velocities from many positions along the line of sight. 

The vertical motions we have discovered might even make a 
planet detectable in a face-on disc, even if it is poorly resolved. In 
this case, the vertical motions induced by the (unresolved) spiral 
arms might well be (mis-)interpreted as turbulent motions. Accord- 
ing to a-disc theory, the turbulent velocity is u tul b ~ a 1,/2 c 3 . If 
the disc were assumed to have a 'normal' value of a w 10 -3 , 
then the sound speed in the disc would be over-estimated by a 
factor of Aiz/ot 1 ^ 2 w 10, implying that the temperature would 
be over-estimated by a factor of roughly three. At this point, we 
would only have a disc with an unreasonably large aspect ratio. 
However, the temperature of the disc might be independently com- 



5 CONCLUSIONS 

In this paper, we have presented the results of 3d numerical ex- 
periments of disc-planet interactions. We have concentrated on the 
case of 1 and 2 Mj planets. Our computations have shown strong 
vertical stirring of disc material by the planet. This stirring has a 
number of potentially observable consequences. 

The gap edge is likely to appear relatively thick for a num- 
ber of reasons. Direct illumination from the star will warm the gap 
edge, slightly increasing the scale height. Density enhancements in 
the spiral arms might also be intepreted as a larger scale height in 
an unresolved image. However, it is the coupling between the gas 
and dust which is most significant. Vertical stirring by the planet's 
spiral arms will counteract dust sedimentation in the disc. Com- 
bined, these effects will make the region close to the planet appear 
"puffed up" relative to the disc. We argue that this strengthens the 
case for a planetary cause of the hole seen in CoKuTau/4, since the 
disc edge is not only sharp, but also quite thick. 

We have also examined the possibility of detecting the ve- 
locity structure produced by the planet using the latest generation 
of telescope, such as ALMA. Rescaling the orbit of our planet to 
10 au, we believe that the presence of a planet is just detectable. 
Even if the disc were unresolved, the effects of a planet might still 
be detectable. A high 'turbulent' temperature in a disc with a cold 
'chemistry' temperature would suggest the presence of a planet. 

There are a number of avenues for further work, which we 
are actively pursuing. Firstly, we are expanding our library of mod- 
els, testing planets of different mass, in discs of varying viscosity 
and aspect ratio. We have demonstrated that planet-induced veloc- 
ity structures should be detectable, but we are not yet in a posi- 
tion to interpret real observations. Compiling a library of models 
will provide us with the tools necessary to interpret future obser- 
vations in detail. Our simulated observations are somewhat crude, 
since they take no account of disc chemistry (which directly af- 
fects which molecules are available for observation; freeze-out of 
CO could be a problem). We will refine our models in the future, 
computing the disc temperature and chemistry in more detail, and 
comparing these to the temperature profile assumed by FLASH. 

We have shown how planets embedded in circumstellar discs 
create significant features in the vertical structure of those discs. 
These structures were not visible in earlier 2d computations. The 
vertical features created are likely to have detectable consequences, 
raising the exciting possibility of catching planets as they form. 
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Figure Al. The effect of the choice of 7 on the Mach number of the vertical 
velocity on planes with z = r H /4 (top) and r H (bottom) for a q = 10~ 3 
planet in the low viscosity disc. The scale ranges from M z = —1 (black) 
to 1 (white). The numerical experiment producing this Figure was identical 
to that of the middle row of Figure|2] except this Figure uses 7 = 1.01 in 
the Riemann solver 
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Figure Bl. The effect of including all components of the viscous stress 
tensor on the Mach number of the vertical velocity on planes with z = f H /4 
(top) and r H (bottom) for a q = 10~ 3 planet in the low viscosity disc. 
The scale ranges from A4 Z = —1 (black) to 1 (white). The numerical 
experiment producing this Figure was identical to that of the middle row of 
Figure[2] except this Figure for the inclusion of the extra viscous terms 



APPENDIX A: EFFECT OF THE EQUATION OF STATE 

All the numerical experiments presented above used an adiabatic 
gas with 7 = 1.1 in the Riemann solver. However, we imposed a 
temperature structure, and did not evolve the energy equation. This 
is not entirely self-consistent. To check the effect of this inconsis- 
tency on our results, we re-ran the low viscosity case for a Jupiter 
mass planet, with 7 = 1.01. 

In Figure lATl we show the effect of this change on the value of 
A4 Z on slices with z — r H /4 (top) and r H (bottom). Apart from the 
value of 7, the numerical experiment used to produce this Figure 
was identical to that of the middle row of Figure [2] We see that 
similar behaviour occurs in the disc with a 7 = 1.01 gas. 

Ideally, one shou ld solve the energy equation in the disc (e.g. 
iD'Angelo etai]|2003l) . perhaps even considering radiative transfer. 
However, this test reassures us that the vertical motions induced 
by the planet are real, and not merely artifacts of the numerical 
method. 



APPENDIX B: ADDING COMPONENTS OF THE 
VISCOUS STRESS TENSOR 

As noted in Section[2] the runs presented above only used the r—cj) 
component of the viscous stress tensor. We made this decision be- 
cause the true nature of viscosity in accretion discs is not known, 
and we simply wished to include a term which would cause angu- 
lar momentum transport. However, our code implements all com- 
ponents of the viscous stress tensor, and in this section, we present 
a run where all components were included. 

In Figure IbTI we show the effect of including all the compo- 
nents of the viscous stress tensor on the value of A4 Z on slices with 
z — r H /4 (top) and r H (bottom). Apart from the inclusion of these 
terms, the numerical experiment used to produce this Figure was 
identical to that of the middle row of Figure [2] We can see that the 
inclusion of the extra viscous terms has not suppressed the strong 
vertical motions we saw previously. As with the case of increased 
viscosity, this is not surprising: the viscosity affects the structure 
of the disc over many orbits, whereas the planet stirs the disc on 
its orbital timescale. Viscosity can act to damp some of the fine 
structure, but it cannot prevent the planet exciting strong vertical 
motions within the disc. 
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Figure CI. The effect of doubling the grid resolution on the Mach number 
of the vertical velocity on planes with z = r H /4 (top) and r H (bottom) 
for a q = 10 -3 planet in the low viscosity disc. The scale ranges from 
Mz = — 1 (black) to 1 (white). This is a continuation of the numerical 
experiment which produced the middle row of Figure[2] with the resolution 
doubled 



APPENDIX C: EFFECT OF RESOLUTION 

As noted in Section[3] our resolution (particularly the vertical reso- 
lution) in the innermost portions of the disc is not especially good. 
This is a result of necessity: global 3d numerical experiments are 
extremely computationally expensive. However, it is important to 
verify that numerical results are not purely artifacts of resolution. 
Accordingly, we instructed FLASH to continue the run with Jupiter 
in a low viscosity disc, but with the resolution doubled. Figure lCTl 
shows the result of this continuation. 

We can see that the structure in M z is retained. However, the 
noise in the inner portions of the disc has been substantially sup- 
pressed (albeit not eradicated), especially in the z — r H plane. The 
peaks in A4 Z appear to be slightly smaller in Figure [Cl] indicat- 
ing that our runs at standard resolution might be slightly under- 
resolving the vertical stirring. This said, the structures seen are very 
similar. 
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